# Goal: Measure effect of social network usage on political participation
# Dependency: data_taps.Rdata

library(MatchIt)
library(weights)

## Load data

load("data_model.Rdata")

#-----------------------------#
#------- Naive ATE  ----------#
#-----------------------------#

 ## Without survey weights:

ATE.naive <- wtd.t.test(data.model$engagement2[data.model$ssn == 1], data.model$engagement2[data.model$ssn == 0], weight = rep(1, sum(data.model$ssn)), weighty = rep(1, sum(1 - data.model$ssn)), alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

 ## With survey weights:

wATE.naive <- wtd.t.test(data.model$engagement2[data.model$ssn == 1], data.model$engagement2[data.model$ssn == 0], weight = data.model$jun2012wt1[data.model$ssn == 1], weighty = data.model$jun2012wt1[data.model$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

#------------------------------------------#
#------- Propensity Score Model  ----------#
#------------------------------------------#

## Without controlling for survey weights

z.out1.ps <- glm(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess,data=data.model, family=binomial(link="logit"))
summary(z.out1.ps)

PS.mle1 <-z.out1.ps$fitted.values

## Controlling for survey weights

z.out2.ps <- glm(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess + jun2012wt1,data=data.model, family=binomial(link="logit"))
summary(z.out2.ps)

PS.mle2 <-z.out2.ps$fitted.values

#--------------------------------------------#
#------- Propensity Score Weighting----------#
#--------------------------------------------#

# Propensity score weights

data.model[data.model$ssn == 1,"PS.weight1"] <- 1
data.model[data.model$ssn == 0,"PS.weight1"] <- PS.mle1[data.model$ssn == 0]/ (1 - PS.mle1[data.model$ssn == 0])

data.model[data.model$ssn == 1,"PS.weight2"] <- 1
data.model[data.model$ssn == 0,"PS.weight2"] <- PS.mle2[data.model$ssn == 0]/ (1 - PS.mle2[data.model$ssn == 0])

# Treatment effects, without survey weights

 ## Without regression adjustment

ATT.weight <- wtd.t.test(data.model$engagement2[data.model$ssn == 1], data.model$engagement2[data.model$ssn == 0], weight = data.model$PS.weight1[data.model$ssn == 1], weighty = data.model$PS.weight1[data.model$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

 ## With regression adjustment

ATT.weight.reg <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = data.model, weights = PS.weight1))$coefficients[2,c(1,2)]

# Treatment effects, with survey weights

data.model[,"adj.PS.weight"] <- data.model$PS.weight2 * data.model$jun2012wt1

 ## Without regression adjustment

ATT.wweight <- wtd.t.test(data.model$engagement2[data.model$ssn == 1], data.model$engagement2[data.model$ssn == 0], weight = data.model$adj.PS.weight[data.model$ssn == 1], weighty = data.model$adj.PS.weight[data.model$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

 ## With regression adjustment

ATT.wweight.reg <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = data.model, weights = adj.PS.weight))$coefficients[2,c(1,2)]

#------------------------------------------#
#------- Nearest Neighbor Matching  -------#
#------------------------------------------#

set.seed(1000)

# Matching

 ## Without survey weights

m.out1<-matchit(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess, distance = "logit", data= data.model, method="nearest", m.order="largest", caliper = 0.1, discard="control", ratio = 2, replace = TRUE)

summary(m.out1)

matched.data1 <- match.data(m.out1)

 ## With survey weights

set.seed(1000)

m.out2<-matchit(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess + jun2012wt1, distance = "logit", data= data.model, method="nearest", m.order="largest", caliper = 0.1, discard="control", ratio = 2, replace = TRUE)

summary(m.out2)

matched.data2 <- rbind(data.model[rownames(m.out2$match.matrix),],
                      data.model[m.out2$match.matrix[,1],],
                      data.model[m.out2$match.matrix[,2],]                               
                      )
matched.data2$W.matched.treated <- rep(data.model[rownames(m.out2$match.matrix),]$jun2012wt1,3)
matched.data2$weights <- match.data(m.out2)[rownames(matched.data2),]$weights

# Treatment effects

 ## Without survey weights

  ### Without regression adjustment

ATT.nn <- wtd.t.test(matched.data1$engagement2[matched.data1$ssn == 1], matched.data1$engagement2[matched.data1$ssn == 0], weight = matched.data1$weights[matched.data1$ssn == 1], weighty = matched.data1$weights[matched.data1$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

  ### With regression adjustment

ATT.nn.reg <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = matched.data1))$coefficients[2,c(1,2)]

 ## With survey weights

  ### Without regression adjustment

matched.data2$weights2 <- matched.data2$W.matched.treated * matched.data2$weights

ATT.wnn <- wtd.t.test(matched.data2$engagement2[matched.data2$ssn == 1], matched.data2$engagement2[matched.data2$ssn == 0], weight = matched.data2$weights2[matched.data2$ssn == 1], weighty = matched.data2$weights2[matched.data2$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

  ### With regression adjustment

ATT.wnn.reg <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = matched.data2, weights = weights2))$coefficients[2,c(1,2)]

#-------------------------------------------#
#------- Subclassification Matching  -------#
#-------------------------------------------#

n.sub <- 8

# Matching

 ## Without survey weights

m1.out<-matchit(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess ,distance = "logit", data = data.model, method="subclass", discard = "control", subclass = n.sub, sub.by = "all")

matched.data1 <- match.data(m1.out)

 ## With survey weights

m2.out<-matchit(ssn ~ dailyiuse + female + agegroup + edushort + incshort + iaccess + jun2012wt1,distance = "logit", data = data.model, method="subclass", discard = "control", subclass = n.sub, sub.by = "all")

matched.data2 <- match.data(m2.out)

# Subclass weights

 ## Without survey weights

subclass.weights1 <- aggregate(matched.data1$ssn[matched.data1$ssn == 1], by =list(matched.data1$subclass[matched.data1$ssn == 1]), sum)$x

 ## Without survey weights

subclass.weights2 <- aggregate(matched.data2$jun2012wt1[matched.data2$ssn == 1], by =list(matched.data2$subclass[matched.data2$ssn == 1]), sum)$x

# Treatment effects

 ## Without survey weights

  ### Without regression adjustment

ATT.s <- wtd.t.test(matched.data1$engagement2[matched.data1$ssn == 1], matched.data1$engagement2[matched.data1$ssn == 0], weight = matched.data1$weights[matched.data1$ssn == 1], weighty = matched.data1$weights[matched.data1$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

  ### With regression adjustment

ATT.s.reg.list <- list()

for (s in 1:n.sub) {
  ATT.s.reg.list[[s]] <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = matched.data1[matched.data1$subclass == s,]))$coefficients[2,c(1,2)]
}

ATT.s.reg <- NULL
ATT.s.reg[1] <- weighted.mean(matrix(unlist(ATT.s.reg.list), ncol = 2, byrow = TRUE)[,1], w = subclass.weights1)
ATT.s.reg[2] <- sqrt(weighted.mean((matrix(unlist(ATT.s.reg.list), ncol = 2, byrow = TRUE)[,2])^2, w = subclass.weights1))

 ## With survey weights

  ### Without regression adjustment

matched.data2$weights2 <- matched.data2$jun2012wt1 * matched.data2$weights

ATT.ws <- wtd.t.test(matched.data2$engagement2[matched.data2$ssn == 1], matched.data2$engagement2[matched.data2$ssn == 0], weight = matched.data2$weights2[matched.data2$ssn == 1], weighty = matched.data2$weights2[matched.data2$ssn == 0], alternative="two.tailed", samedata = FALSE)$additional[c(1,4)]

  ### With regression adjustment

ATT.ws.reg.list <- list()

for (s in 1:n.sub) {
  ATT.ws.reg.list[[s]] <- summary(lm(engagement2 ~ ssn + dailyiuse + female + agegroup + edushort + incshort + iaccess, data = matched.data2[matched.data2$subclass == s,], weights = weights2))$coefficients[2,c(1,2)]
}

ATT.ws.reg <- NULL
ATT.ws.reg[1] <- weighted.mean(matrix(unlist(ATT.ws.reg.list), ncol = 2, byrow = TRUE)[,1], w = subclass.weights2)
ATT.ws.reg[2] <- sqrt(weighted.mean((matrix(unlist(ATT.ws.reg.list), ncol = 2, byrow = TRUE)[,2])^2, w = subclass.weights2))

# Table 2

rbind(ATE.naive, wATE.naive, ATT.weight, ATT.wweight, ATT.weight.reg, ATT.wweight.reg, ATT.nn, ATT.wnn, ATT.nn.reg, ATT.wnn.reg, ATT.s, ATT.ws, ATT.s.reg, ATT.ws.reg)


